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Study of gravitational radiation from cosmic domain walls 
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In this paper, following the previous study, we evaluate the spectrum of gravitational wave back- 
ground generated by domain walls which are produced if some discrete symmetry is spontaneously 
broken in the early universe. We apply two methods to calculate the gravitational wave spectrum: 
One is to calculate the gravitational wave spectrum directly from numerical simulations, and an- 
other is to calculate it indirectly by estimating the unequal time anisotropic stress power spectrum 
of the scalar field. Both analysises indicate that the slope of the spectrum changes at two charac- 
teristic frequencies corresponding to the Hubble radius at the decay of domain walls and the width 
of domain walls, and that the spectrum between these two characteristic frequencies becomes flat or 
slightly red tilted. The second method enables us to evaluate the GW spectrum for the frequencies 
which cannot be resolved in the finite box lattice simulations, but relies on the assumptions for the 
unequal time correlations of the source. 

PACS numbers: 98.80.Cq, 04.30.Db 



I. INTRODUCTION 



Gravitational wave (GW) is one of the robust predic- 
tions of general relativity, and expected to be detected 
in the next decades. Since GWs have few interactions 
with matter and radiation, they propagate almost freely 
after their production. Therefore, analogously to the cos- 
mic microwave background, the search for a stochastic 
background of GWs will give us rich informations about 
the early universe which has not been probed by electro- 
magnetic waves. Various mechanisms to generate GW 
backgrounds have been proposed, such as inflation 1|, 
(p)reheating after inflation 0-111) cosmic strings 0, Hj, 
and first order phase transitions 0-|l3- addition to 
them, we propose that domain walls, which are surface 
like topological defects produced when a discrete sym- 
metry is spontaneously broken, can be another source 
of the stochastic background of GWs. The existence of 
domain walls is cosmologically unacceptable, since they 
eventually overdose the universe [ll]. However, if do- 
main walls are unstable and decay at a sufficiently early 
time [3, [l5| , the energy stored in them would be radi- 
ated as GWs, which become the stochastic background 
observed today [il,|T3| . There are several particle physics 
models which predict such phenomena. For instance, the 
spontaneous breaking of the discrete R symmetry in the 
theory with supersymmetry induces domain walls which 
decay when the Hubble parameter becomes comparable 
to the scale of the gravitino mass |18i, ilQ]. Also, the the- 
ory of axions, which was introduced in order to solve the 
strong CP problem of quantum chromodynamics, nat- 
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urally predicts the existence of domain walls [20|, |2lj. 
Therefore, the observation of the stochastic background 
of GWs produced by domain walls can become another 
probe of the theory beyond the standard model of parti- 
cle physics. 

The stochastic background of GWs of cosmological ori- 
gin is expected to be isotropic, stationary and unpolar- 
ized, and therefore characterized by its frequency spec- 
trum In the previous study [Ol, we calculated the 
spectrum of GWs produced by domain walls based on the 
numerical simulation of the scalar field in the expanding 
universe. The results were straightforwardly obtained, 
and it was shown that GWs from domain walls have a 
broad and nearly flat spectrum. However, in the pre- 
vious work, we chose the energy scale of the symmetry 
breaking as an unrealistic value rj ~ 10"'^''GeV in order 
to follow the domain wall evolution from the initial ther- 
mal state, and simply extrapolated the numerical result 
to predict the spectrum observed today. This estima- 
tion gives the large uncertainty, which can be a factor 
of 0(10^1) in the magnitude of GWs. Furthermore, we 
put the thermal initial condition for numerical simula- 
tions, which made an additional peak at high frequencies 
in the GW spectrum. This might contaminate the spec- 
trum of GWs purely generated from domain walls. In 
order to remove these difficulties and give more accurate 
predictions for future observations, it is necessary to in- 
vestigate further about the spectrum of GWs. 

In this work, we apply two approach to evaluate the 
spectrum of GWs produced by domain walls. First, we 
perform three dimensional lattice simulation of domain 
walls, and calculate the GW spectrum directly from the 
results of the numerical simulations. This procedure is 
similar to that in our previous work, but we work with 
different setups: While we performed the calculations 
only for radiation dominated background in the previous 
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work, we calculate the GW spectrum both in radiation 
and matter dominated backgrounds. Furthermore, we 
do not assume the coupling with the thermal bath. In 
this case, the GW spectrum would represent the feature 
of that purely produced by domain walls. This analysis 
will confirm and clarify the features of the GW spectrum 
which we found in the previous work. Second, we reevalu- 
ate the GW spectrum indirectly by using the anisotropic 
stress power spectrum obtained from numerical simula- 
tions. For the evaluation of the GW spectrum, we intro- 
duce the same approximations which has been used to 
calculate the spectrum of GWs generated by first order 
phase transitions [23l - [25| . This method depends on vari- 
ous assumptions about the source of GWs, but it might 
become an alternative way which enables us to evaluate 
the GW spectrum without relying on numerical simula- 
tions with a long dynamical range. 

This paper is organized as follows. In section |lll we 
describe a notation to calculate GWs and derive a for- 
mula which gives the spectrum of GWs produced in the 
arbitrarily expanding background. Then, we present the 
results of numerical simulations and discuss the features 
of the spectrum of GWs produced by domain walls in 
section Hill In section ITVl we take another approach and 
reevaluate the GW spectrum. We also give forecasts for 
future observations based on the result of the analysis 
performed there. Finally we conclude in section IVl 



II. A STOCHASTIC GRAVITATIONAL WAVE 
BACKGROUND IN THE EXPANDING 
UNIVERSE 

In this section, we derive basic equations for GWs pro- 
duced in arbitrary expanding background. We consider 
a spatially flat Friedmann- Robertson- Walker background 
in which GWs are represented by the spatial metric per- 
turbation 



(1) 



where hij satisfies the transverse-traceless (TT) condi- 



tion dihi 



0. Since we will investigate the gen- 



eration of GWs in both radiation and matter dominated 
background, it is convenient to consider the arbitrarily 



expanding background where the scale factor evolves as 

a(t) cx r" cx i'^, (2) 

where t is confornial time defined by dr — dt/a. 

The metric perturbations hij obey the linearized Ein- 
stein equation 



hij (t , x) -I- 3Hhij {t, x) 



-hij{t,ic) 



167rG, 



(3) 



where a dot denotes a derivative with respect to cosmic 
time t and T'^^ is the TT part of the stress-energy ten- 
sor. If we work in spatial Fourier space and change time 
variable from cosmic time t into conformal time r, this 
equation gives 

h'lj (r, x) + ^ /i^ (r, k) -t- fc2 /i^^. (r, k) 

where a prime denotes a derivative with respect to con- 
formal time T. Defining the rescaled metric 



ah 



(5) 



we obtain 



hij{T, k) 



167rG'a(T) TT 



where x = kr, and i' is defined by 

1 3/3-1 



2(1-/3) 



(6) 



(7) 



Let us assume that the source term T,^""" is nonzero 
during the interval < t < Tf. The solution of eq. (joj) 
with initial conditions hij{Ti) = h[j{Ti) = is given by 
the time integral of the source term convoluted with a 
Green function 



J 



/ly (T, k) 



dy{yxY/^[N,{x)My) - Mx)N,{y)]a{y)Tj;^ [y,]^) (for r < Tf) 



(8) 



where (x) and [x) are Bessel function and Neumann 
function, respectively. Note that this is just a generaliza- 
tion of the Green function solution obtained in 5] for the 
radiation dominated background. Substituting v — 1/2, 



one can easily check that eq. (|8]) reduces to the result in 
the radiation dominated universe derived in Q . 

After the time Tf, the source term becomes negligible 
in eq. © , and hij is given by a linear combination of two 
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independent solutions of eq. ([6]) without the source term 
/iy(r,k) = A,j(k)(fcr)i/2j,(fcT) 

+ By(k)(fcT)l/2iV,(fcT) {i0VT>Tf). (9) 

The coefficients Aij and Bij are determined by matching 
the solution given by eq. ([8]) with eq. ([9]) a.t t — Tf. We 
obtain 

Ay(k) = / ' dxy^aix)N,{x)T^^^{x,k), 



fc2 

By(k) = -p-y dTv^a(x)J^(x)TTT(a;^k). (10) 

Let us define a dimensionless anisotropic stress tensor 

«"'WTTT(r,k) = (p + p)%(r,k), (11) 

where p is the background homogeneous energy density 
and p is the background homogeneous pressure. We as- 
sume that the source is statistically homogeneous and 
isotropic, and introduce the unequal time correlator of 
the anisotropic stress tensor 

(%(Ti,k)n^(T2,k')) = (2^)3(5(3) (k-k')n(fc,ri,T2), 

(12) 

where (• • • ) denotes an ensemble average for a stochastic 
background. We note the following relations 



THa{T) — a 



3(1 + w)' 



(13) 
(14) 



where ui is a mean equation of state defined hy p = wp. 
By using these relations, we can rewrite eq. ()10p as 



i?,,(k) 



27r/3 
(W)2 J.. 



-J 

dxx~^''^a{x)Nt,{x)Ilij{x,k), 
dxx^^^^a{x)J^{x)Ilij {x, k). 



(15) 



The energy density of GWs is given by (see e.g. ^]) 



1 



32^G ' 
1 

327rGa4(7 



■{KAt, x)/i^ -(t,x)), 



(16) 



where we neglected the terms with higher order in aH in 
the second equality, since we assume that the wavelength 
of GWs is well inside the Hubble radius at the time r 
(fcr » 1). Substituting eqs. (O and (IT5|) into eq. pB]) . 
and using eq. p^. we obtain 



Pgw{t) 



327rGa4(t) J (27r)3 tt (1 - /S)-^ 

^/ dxi dx2 

a(xi)N,y{xi) / ^r-:-a(x2)iVi,(a;2)n(fc, ri,T2) 



Xi '- 



where we used the approximations for /cr 3> 1 



/ -575-a(a;2)J^(a;2)n(A;,Ti,T2) > 
Jxi x^ J 



a{xx)Ji,[xx) \ -^a{x2)Jv{x2)^i}l,Tx,T2)\ , (17) 



J.(fcr) ^ cos (fcr - ^ - J) , iV.(fcr) ^ -n (fcr - ^ - J) , (18) 



r 



and averaged over a period of the oscillation of sine and density of GWs at the time t as 
cosine with time. We define the fraction of the energy 

1 dpg„{t ) 
p(t) d\nk 



^.At) = -zk-^. (19) 
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Note that, it is not the value which would be observed at 
the present time. We introduce this notation for conve- 
nience to present the result of numerical simulations. We 



J 



will convert it into the spectrum of GWs observed today 
in section ITVDI Substituting eq. (IT71) into eq. (ITOl) . we 
finally obtain 



fc3 



67r(l - /3)2 



dxi X a{xi) 
'xi xi a{x) 



dx2 X a(x2) 
2^X2 a{x) 

I 



[N,{XI)N,{X2) + Mx,)J,{x2)]U{k.Tl,T2), (20) 



where we used p{t) — ■^^H'^{t). 

In the previous work the spectrum of GWs was 
directly obtained by calculating TT pert of the stress- 
energy tensor in numerical simulations. Another 
way to evaluate the spectrum of GWs is to estimate the 
anisotropic stress power spectrum n(fc,Ti,T2) and use 
eq. (I^H)) . We will evaluate in both ways and compare 
the results. The advantage of using eq. ((20|) is that wc 



this effect by adding a term 



can decompose the origin of k dependence of figw: The 
tilt of the spectrum is determined by the k dependence in 
n(fc, ri,T2) and the time integral of the Green function. 
This decomposition might be helpful to understand the 
precise form of the GW spectrum, as we see in section lTVl 



5V = erj(j) ( -0^ — if 



(23) 



to the potential given by eq. (|22p . The dimensionless 
parameter e, which we call "bias" , controls the magnitude 
of energy difference between two vacua and determines 
the life time of domain walls. We assume that e is much 
smaller than 1 since we are interested in the circumstance 
in which the discrete symmetry is held approximately. In 
particular, the condition e < 0.15A must be satisfied in 
order that the infinite size of domain is formed [l7j . 



III. NUMERICAL SIMULATION OF DOMAIN 
WALLS 

In this section, we show the results of numerical simu- 
lations. We consider the simple model of real scalar field 
(j) in which a discrete Z2 symmetry is spontaneously bro- 
ken. The evolution of in the expanding background is 
described by the Klein-Gordon equation 



3iJ0- 



where the potential is given by 



A, 



dV 



0, 



2\2 



(21) 



(22) 



This potential has Z2 symmetry under which the scalar 
field transforms as </> — > —4>. This symmetry is sponta- 
neously broken and scalar field gets vacuum expectation 
value {(p) = ±rj. Domain walls can be formed around the 
region where the value of the classical field changes from 
— r/ to +ri. 

If such domain walls were created in the early uni- 
verse and survived until today, they eventually come to 
overdose the energy density of the universe and disturb 
the success of standard cosmology [ll] . One way to avoid 
this problem is to introduce a term in the potential which 
explicitly breaks the discrete symmetry and lifts the de- 
generacy of vacua If such a term exists, walls 
become unstable and eventually disappear. We model 



A. Initial conditions 

Eq. ([2T|) . which describes the evolution of domain 
walls, is highly nonlinear and difficult to solve analyti- 
cally. Then we solve eq. (PT|) numerically on the three di- 
mensional lattice. In the previous study 



17| , scalar field 



is assumed to be in thermal equilibrium with tempera- 
ture T, and the initial field configurations are generated 
by considering finite temperature effects. In this primary 
thermal stage, the scalar field fiuctuations also produce 
GWs with the spectrum peaked at the frequency corre- 
sponding to the mass of the scalar field ~ ^fXrj. However, 
this setup might be irrelevant to our interest to calculate 
GW spectrum produced by domain walls. The reason is 
as follows: The amplitude of GWs becomes large enough 
to observe if domain walls survived for sufficiently long 
time. This means that the spectrum of GWs produced 
at the primary stage is negligible compared with that 
produced by domain walls at the late time. Further- 
more, the assumption that the phase transition occurred 
at T ~ ?7 gives a severe constraint on the range of param- 
eters which we choose to perform realistic simulations. In 
numerical simulations, we must resolve the width of the 
wall (5iL, ^ rj^^ and keep Hubble radius smaller than 
the size of the simulation box. If we assume that the tem- 
perature is given by T ~ 7y, the ratio of these two length 
scale is Sw/H~^ ~ rj/Mp, where Alp is the Planck mass. 
Therefore, we must choose rj close to the Planck scale in 
order to maintain the resolution of the width of domain 
walls. In numerical simulations performed in [l7| . the 
value of 1] is chosen to be 10^''GeV, and the results of 
the simulations are extrapolated into lower values of rj. 
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However, it is not obvious that this extrapolation is held 
for arbitrary scale of 77. 

To avoid these difficulties, we omit the assumption of 
thermal initial conditions and give the initial field con- 
figurations as Gaussian random amplitudes. In this case, 
we expect that only domain walls contribute as a source 
of GWs. In addition, if we normalize all the dimensionful 
quantities in the unit of 77, the results of numerical sim- 
ulations become independent of rj. See appendix for 
more details of the setup of the simulations. 



B. Evolution of the domain wall networks 

Now, we present the results of numerical simulations. 
We performed lattice simulations with 256^ points in 
both radiation and matter dominated backgrounds. The 
comoving size of the simulation box is set to be 50 (in 
the radiation dominated era) or 24 (in the matter dom- 
inated era) in the unit of 7/^^. We fix A = 0.1 and vary 
the value of e. The initial time (in the unit of 7y~^) is 
set to be — 1. The final time is set to be i/ — 151 
in the simulation with radiation dominated background 
and tf = 65 in the simulation with matter dominated 
background. 

Figure [1] shows the time evolution of the area occupied 
by domain walls in the simulation box. We see that, if 
e 7^ 0, the area density of domain walls decays at late 
time. It means that domain walls collapse due to the 
existence of the bias. On the other hand, if e = 0, the 
comoving area density of domain walls evolves as cx t~^. 
This property is called the scaling solution (27l - [3l| and 
corresponds to the fact that the energy density of domain 
walls evolves like pwaii ^ cr/t, where a = 2V2Xr]'^/3 is the 
surface mass density of the domain wall. Figure shows 
the time evolution of the energy density of the scalar 
field. From this figure, we see that the scaling regime in 
which energy densities evolve as oc 1/t begins at around 
t ~ 20. 

Note that, in our study, the dynamical range of the 
simulation is quite short. This is due to the fact that 
we can not resolve the width of domain walls for a long 
time in the comoving box if we take account of the cos- 
mic expansion. In particular, the actual dynamical range 
in the simulation with matter dominated background is 
as smaU as (i//tform)"^^^ — (65/20)^/^ ~ 1.8 in conformal 
time, where iform is the time when domain walls enter 
the scaling regime. Future simulations with higher spa- 
tial resolutions should improve the dynamical range and 
confirm our current results of numerical simulations. 



C. Spectrum of gravitational waves 

We calculate the spectrum of GWs directly from nu- 
merical simulations. The method of the calculation is 
summarized in Appendix |B] Here we present the results 
of the direct calculations and briefiy discuss the features 



of the spectrum of GWs. We will reevaluate the spectrum 
of GWs in the next section by using another formalism. 

The spectra of GWs obtained from simulations are 
shown in figure [31 The vertical axis represents the am- 
plitude of GWs defined by eq. ([19]). We normalize it by 
the dimensionless quantity 



piU) 3(3- 



(24) 



where = Grj^ is the energy density of GWs (esti- 
mated by using the quadrupole formula of GWs) radi- 
ated by a source which has a characteristic scale i], and 
p{U) = 3H^{U)/87rG = 3^2^V87rG is the background 
homogeneous energy density at the initial time of the sim- 
ulations. We introduce this notation for a convenience to 
present the numerical results: r2gw(i)/^i; becomes C(l) 
at the beginning of the simulation. The horizontal axis 
represents the comoving wavenumber k normalized by 
rj. This is not the frequency of GWs, which is given by 
/ k/2TTa{t). 

From figure [31 we see that the spectrum of GWs is al- 
most fiat in both radiation dominated and matter domi- 
nated backgrounds. In the previous study it is con- 
jectured that this fiat spectrum extends roughly between 
the frequency corresponding to the Hubble radius at the 
time of the decay of domain walls and that correspond- 
ing to the width of domain walls. Let us define comoving 
wavenumbers kh and fc^, which correspond to the Hubble 
radius and the width of domain walls, respectively. 



^ = 2nH{t) 
a{t) 



and 



a{t) 



(25) 



Note that the values of fc/j and fc„ vary with time. In the 
present numerical simulations where all of the dimension- 
ful quantities are normalized in the unit of 77, kh changes 
from TT at the initial time to O.OStt at the final time in the 
simulation with radiation dominated background, and it 
changes from I.Stt to O.SStt in the simulation with matter 
dominated background. Also, k^ changes from O.GStt to 
T.Stt in the simulation with radiation dominated back- 
ground, and it changes from 0.637r to IOtt in the sim- 
ulation with matter dominated background. Therefore, 
these two scales become separate at the late time of sim- 
ulations. In fact, we see that the band of the spectra 
shown in figure [31 becomes wider as time passes. 

In our previous study [l^l, the spectrum of GWs had 
a peak at the high frequency (around k < kyj), and the 
slope of the spectrum was mildly blue tilted. However, 
in the present results, the spectrum with e = becomes 
slightly red tilted as we see in figures [31 (a) and (c). This 
difference might be caused by the fact that we use dif- 
ferent initial conditions for numerical simulations: The 
spectrum obtained in [17,] contains the GWs produced 
during the primary thermal stage (see section IIII Ap . 
which is not completely negligible in the time scale of 
the simulations and gives a peak at the high frequency. 
Therefore, we expect that if GWs from primary ther- 
mal stage become negligible and the source of GWs is 
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FIG. 1: The time evolution of the comoving area density of domain walls for various values of e in (a) radiation dominated 
background and (b) matter dominated background. 





FIG. 2: The time evolution of the kinetic/gradient/potential energy densities of domain walls (in the unit of 77*) in (a) radiatL 
dominated background and (b) matter dominated background. 



purely domain walls, the spectrum of GWs is almost flat 
and slightly red tilted. Then, if domain walls collapse, 
the amplitude of GWs grows at the high frequency as 
we see in figures [3] (b) and (d). This corresponds the fact 
that the false vacuum regions become small when domain 
walls decay. 



To summarize, the spectrum of GWs obtained by di- 
rect numerical calculations indicates that the shape of 
the spectrum is nearly flat (or slightly red tilted) and 
the slope of the spectrum seems to change at two char- 
acteristic frequencies corresponding to the Hubble radius 
at the decay of domain walls and the width of domain 
walls. High frequency modes grow at the late time due 
to the collapse of domain walls. These properties are 
seen in both radiation dominated and matter dominated 
backgrounds. 



IV. THE SPECTRUM OF GRAVITATIONAL 
WAVES FROM DOMAIN WALLS: ANOTHER 
APPROACH 



The lattice simulations revealed that some properties 
of GWs produced by domain walls, such as the nearly 
flat spectrum and the existence of two relevant scales 
(the Hubble radius at the decay of domain walls and 
the width of domain walls). However, these results have 
many numerical uncertainties especially in the high and 
low frequency bands of the spectrum, because of the lack 
of the spatial resolution in the lattice simulation. In order 
to remove this difficulty, we take another approach to 
evaluate the spectrum of GWs and discuss further about 
the shape of the spectrum of GWs. 

In section [ill '^e derived eq. ([20| . which indicates 
that we can evaluate the GW spectrum if we know the 
anisotropic stress power spectrum H(/c,Ti,r2). It is diffi- 
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FIG. 3: The spectra of gravitational waves obtained by numerical simulations. The top panel shows the results in radiation 
dominated background with (a) e = and (b) e = 0.003, and the bottom panel shows that in matter dominated background 
with (c) e = and (d) e = 0.004. The different colors correspond to the spectra at different time. The spectra are shown from 
the time t — 31 (pink) tot= 151 (green) with the interval At — 20 in the simulation with radiation dominated background, and 
from the time t = 25 (pink) to t = 65 (green) with the interval At = 8 in the simulation with matter dominated background. 



cult to evaluate unequal time power spectrum n(fc, ti, T2) 
without relying some ansatzs about the unequal time cor- 
relation of the source of GWs. In section IIV Al we de- 
scribe approximations to evaluate the unequal time corre- 
lation functions, which is used to estimate the GW spec- 
trum from first order phase transitions in the literature. 
Using these approximations, we reevaluate the GW spec- 
tra and compare them with the results of direct numerical 
calculations in section HV Bl Then, we comment on the 
property of the GW spectrum obtained from these analy- 
sises in section HV CI Finally, in section HVDl we convert 
the results into the present-day observables. 



Approximations for the unequal time 
correlation functions 



this quantity is given as a product of T^j (t, k) , we hope 
that it also be computed directly from lattice simulations. 
To be exact, we can not compute the ensemble average 
in the left hand side of eq. (|12l) from a single realization 
of the numerical simulation. However, we approximately 
calculate the ensemble average by taking average over a 
large volume. Thus, 

J ^(n,, (k,r)n^(k,r))« J ^n„ (k,r)n^(k,r) 

(26) 

is satisfied in the limit V — > 00, where Jdh/Air is an 
average over the directions of k, and V is the comoving 
volume of the simulation box. It is easy to show eq. p6|) , 
assuming Gaussian statistics (see e.g. (32j). Using the 
approximation given by eq. (j26p . we obtain 



We would like to evaluate n(/c,ri,T2) defined by 
eq. (1121) . Before estimating this unequal time correlator, 
let us consider the equal time correlation function. Since 



n(/c,r, r) 



167r2G2/?2 fd^k 



a^{T)H^V J 47r 



-Tr{T,\.)T^*{T,\.) 



(27) 
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where we replaced the factor {2TrfS'-^^0) ~ V. By using 
eq. (P7)) . we can compute the equal time power spectrum 
n(fc, r, r) directly from lattice simulations. 

Next, consider the unequal time correlator n(fc, ti, T2). 
In the recent study of GW generation from bubble colli- 
sions [23], it was suggested to use some approximations 
for the unequal time correlator in order to evaluate the 
GW power spectrum. These approximations are also 
applied to evaluate GWs from turbulence and inagnetic 
fields generated by a first-order phase transition [2^ (see 
also |24|). 

There are three kinds of approximations: 

1. Totally coherent approximation 

The source at different time is perfectly correlated, 
and the unequal time correlator is given by 

n(fc,ri,r2) - Vn(A:,ri,Ti)Vn(fc,T2,r2). (28) 

2. Incoherent approximation 

The source at different time is not correlated, and 
the unequal time correlator is given by 

n(fc,Ti,r2) =n(fc,Ti,ri)^(ri-T2)Ar, (29) 

where At is a characteristic time scale for the du- 
ration of the source (left as a free parameter). 

3. Top hat approximation 

The source is correlated for modes with a time sep- 
aration |ri — T2I < Xc/k, and the unequal time cor- 
relator is given by 

n(fc, Tl,T2) 

= n(fc,Ti, ri)e(r2 - Ti)e - - n)) 
+ n(fc, T2, T2)e(ri - T2)e - (n - T2)) , (30) 



where dimensionless parameter of 0(1) and 

8(t) is the Heaviside function. 

The first two approximations are physically less moti- 
vated and just introduced for a comparison. On the other 
hand, the top hat ansatz is an intermediate case between 
two extreme cases. The physical interpretation of this 
approximation is that the correlation is lost for a time 
difference larger than about one wavelength. 

Using these approximations, we can evaluate the un- 
equal time correlator n(fc, ti,T2) from the equal time cor- 
relator Tl{k, T, t) obtained by the numerical simulations. 
In the next subsection, we evaluate the GW spectrum by 
combining approximations given by eqs. (|28 |) -p0 |) and 
the formula for the amplitude of GWs given by eq. (|20l) . 



B. Evaluation of the gravitational wave spectrum 

First, let us evaluate the equal time anisotropic stress 
power spectrum n(A;,r, r). By using eq. (j27p . we can 
compute n(fc,r, r) from lattice simulations. The result 
is shown in figure |4l From this figure, we see that the 
tilt of the spectrum becomes steeper in the small scale. 
In section IIII CI we argued that the spectrum of GWs is 
determined by two characteristic scales: the width of the 
wall and the Hubble radius. Therefore, we expect that 
the tilt of the power spectrum Il{k, r, r) also changes at 
these characteristic scales. Furthermore, since the source 
has no correlation for the length scale beyond the Hubble 
radius, the power spectrum would become independent of 
k in the large scale limit (this is just an assumption, but 
one can show analytically that H(A; — > 0) is independent 
of k for a spherical bubble configuration [1^). Based on 
these considerations, we express the k dependence of the 
power spectrum obtained by the numerical simulation at 
the final time Tf as 



U{k,Tf,Tf)/G'v^A 



KUrf) 



D 



(31) 



where Kh{Tf) and Kw{Tf) are given by 



KhiTf) = k/ku{Tf) = 

Kw{Tf) = k/k^{Tf) = 




These are the ratio between the comoving momentum 
k and k^, or fc^,, defined by eq. ()25p at the conformal 
time Tf. We fit the result of H(A:,ry,Tj) obtained from 
numerical simulations to the function (|31l) by using the 



least squares method and determine unknown parameters 
A, B, C, D and E (see appendix [Cl for details). The 
results are shown in table HI The parameters are fixed 
with uncertainties of 0(1) - O(0.1)%. 

We note that the amplitude of n(A:,r, r) grows with 
time, as we see in figure HI To include this property, we 
model the form of the equal time anisotropic stress power 
spectrum as 

U{k,T,T)(x(^^^^ L^Sik,T), (33) 
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TABLE I: The values of the parameters in eq. (|31|l determined by the least squares fitting. A is the amplitude of n(fc, r, T)/G'^'ri 
at T = rj~^, as defined in eqs. p4)l and (f35)) . 





radiation dominated 


matter dominated 




(1.84±0.09)xf0'-' 


(f.57±0.03)xf0'' 


B 


0.388±0.008 


0.724±0.004 


C 


2.020±0.005 


2.583±0.002 


D 


0.f866±0.0007 


0.28f3±0.0008 


E 


3.f448±0.0004 


3.5444±0.0009 


A 


(3.4±0.2)xf0"^ 


(2.9±0.f)xf0"'* 



where pgrad is the gradient energy density of the source, 
L is a characteristic scale of the problem, and 5(fc,T) 
is a dimensionless function of k (and r) ^25,]. Since the 
characteristic scale of domain walls is given by the Hub- 
ble radius, it is natural to expect that L is as much as 
r. The scaling solution also implies that pgrad evolves as 
oc cx r"~^, while p evolves as oc cx t^^ in radiation 
dominated background. Similarly, we expect Pgrad 
and p oc r"^ in matter dominated background. There- 
fore, we expect that the amplitude of the power spec- 
trum evolves as (pgrad/ p)^L^ oc t'^, where 7 = 7 in 
radiation dominated background, and 7 = 9 in matter 
dominated background. Combining it with eq. (j31[) . we 
obtain the following expression for the equal time corre- 
lator n(fc, T, r) 

n{k,T,T)^G'7l{TTjyS{k,T), (34) 

(35) 




where A is a rescaled parameter whose value is given 
m table m Since r is normalized by ^ in numerical 
simulations, we write the r-dependent factor as [tt])^ 
in eq. We also assume that S{k,T) depends on 

r through the functions Kh{T) and iir^(r), which are 
obtained by extrapolating the ratios ([5^ into general 
time T. 

We plot the function ([M)) in figure [S] Comparing fig- 
ures [4] and [5l we see that the function given by eqs. (p4)) - 
([571 indeed reproduces the form of n(A:,r, r) obtained 
from numerical simulations except that the spectrum in 
figure [S] does not agree with that in figure H] at early 
times when domain walls do not enter the scaling regime 
in the numerical simulation. Fortunately, this disagree- 
ment does not significantly affect the final form of the 
GW spectrum because the value of n(fc, r, r) at early 
times is much less than that at the final time by a factor 
of 0(10-4). 

Now we have all the ingredients to evaluate the spec- 
trum of GWs. The amplitude of GWs is given by eq. (^01) . 
Setting /3 = 1/2 and = 1/2 (in radiation dominated 
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k [in unit of Til 
(a) 




k [in unit of n] 



FIG. 5: The expression for n(fc,r, r) given by eqs. p4[) - (|37l l in (a) radiation dominated background and (b) matter dominated 
background. 



background), we obtain 

^ 3 dxi r' 



= ^^fc / / cos(xi - a:2)n(fc, Ti,T2), 

(38) 



where the subscript * represents the fact that it is not 
the spectrum of GWs at the present time but that in 



J 



the radiation dominated era (since the energy density 
of GWs is diluted as pgw oc a~^, the ratio (figw)* — 
Pgw/p is independent of t in the radiation dominated era). 
Combining it with approximations given by eqs. (|28|) - 
(PO)) . we can express the spectrum of GWs as a time 
integration of the function S{k,x) given by eq. (I35p 



1 ^kv 1 



1 /kVAx rf 



87r3 Vryy xl 



dxx^S{k, x) 



1 fky 2 ^ r dy , , 

~ ~^ \ <J-xx° —cos{x-y)S{k,x) 
^ STT-i \r]J xjj J^^ y 



r 



dxx ' cos X ^S{k,x)\ + { / dxx"^/"^ smxy/S{k, x) ) > totally coherent 



incoherent 
top hat, 



m 



where x,j = krj ^, Ax = fcAr, x = min{x/, Xc + x}, and Similarly, setting (3 — 2/3 and = 3/2 in matter dom- 

flfj is given by eq. ([24]). inated background, we obtain 



^(t) ^|(r/^^^'^3/2(a;)v/^(fc:^) +[j^f dxx''jy2{x)yW^: 



totally coherent 
incoherent 



J^f dxxlO [iN,/2ix)r + (J3/2(X))2] ^(fc^^) 

^ /;/ /; dyy^/^[Ny2{x)N3^2iy) + J3/2{x)J3/2{y)]Sik,x) top hat. 



(40) 



Note that, in the matter dominated era, figw in eq. ([20)1 means that eq. PO]) represents the GW spectrum at the 
depends on time. In eq. ([40]) . we fixed r = r/, which 
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conformal time Tf. This is the same quantity which 
we computed directly from lattice simulations in sec- 
tion IIII CI (see appendix [B] for details) . 

In figure El we show the spectrum of GWs evaluated 
from eqs. (p9| and (|40| for each approximation, and also 
the spectrum obtained by direct numerical calculations 
which we described in section llll CI In this figure we fixed 
the values At = 1 and Xc = 1. We found that the spec- 
trum with top hat approximation has an agreement with 
the result directly obtained from numerical simulations 
within a factor of C(l). 



C. Contribution from the time evolution of the 
source function 

Now, let us examine the properties of approximations 
which we used in the previous subsections and consider 
their implications for the shape of the GW spectrum. 

Inspection of eqs. ([M)) and (1^0]) tells us that the k 
dependence of ilgw is determined by three factors: the 

I 



phase space volume k^, the time integral of the Green 
function and the function S{k, r) in the equal time 
anisotropic stress power spectrum. In order to clarify 
the role of the individual terms, let us simply ignore 
the time dependence of S{k, r) and write the equal time 
anisotropic stress power spectrum as 



Il{k,T,T)^G'vf{T)S{k), 



(41) 



where /(r) is a dimensionless function of r. In the pre- 
vious subsection, we used /(t) = (t?^)^ for domain walls, 
where 7 = 7 in radiation dominated background, and 
7 = 9 in matter dominated background. Substituting 
eq. (HI]) and approximations (P5 1) - ([5n|) into eq. (PH)) . we 
obtain the following expression 



1 fk\^ 

ngw(T/)/a, = ^ ( - j F{k,n,Tf 



)S{k), (42) 



where F(k, Ti, Tf) is given by 



F{k,Ti,Tf) 



-Irr. 1-/3 



2(1-/3) 



X < 



1^ _3a_T 



Ax I dxx^-^+^ 



dxx 1- 



|3 _ 3 



dxx 1- 



'^Jy{x) 



{N,{x)f + {J,{x)f 



totally coherent 
incoherent 



(43) 



+7 / dyy^-'i [N,{x)N,{y) + J,{x)J,{y)] top hat. 

I 



In figure [7l we plot F{k,Ti,Tf ) as a function of k for 
each approximation. This result shows that the factor 
F{k,Ti,Tf) gives another contribution to the shape of 
the GW spectrum: The factor F{k,Ti,Tf) is suppressed 
for high frequencies k > t^^ in the totally coherent ap- 
proximation and the top hat approximation, while it is 
independent of k in the incoherent approximation [25| . 
This behavior is understood as a interference of the func- 
tions Nu[x) and J^ix), which rapidly oscillate for the 
frequency larger than t^^. We find that F{k,Ti,T-f) de- 
cays like k^'^ for the coherent case and k~^ for the top 
hat case. This means that in the totally coherent case 
the interference is stronger than that in the top hat case. 
Therefore, there are much suppression at high frequen- 
cies for the totally coherent case, as we see in figure |6l 
On the other hand, the behavior of F{k,Ti,Tf) at low 
frequencies k < tJ^ becomes different between the re- 
sult with radiation dominated background and that with 
matter dominated background. This is understood as 
a difference in the behavior of the Neumann function 
in the limit x — >■ 0. Noting that Ny{x) oc x~'' for 



x <C 1 and the integral in eq. is dominated by the 
contribution around x ~ Xi, one can easily show that 
F{k, Ti, Tf) cx fc^"^" for the modes with k <C tJ"^ . There- 
fore, in radiation dominated background {u — 1/2), the 
amplitude of F{k,Ti,Tf) does not depend on k at lower 
frequencies. This amplitude roughly scales as oc tJ in 
the totally coherent approximation and the top hat ap- 
proximation, and as cx tJ ^At in the incoherent approx- 
imation. This implies that the result with the incoherent 
approximation underestimates the GW amplitude than 
other two by a factor of Ar/r/, due to the presence of 
the additional factor At. By contrast, in matter domi- 
nated background (^ — 3/2), F{k,Ti,Tf) behaves as k^^ 
at low frequencies k < tJ^ . This makes the tilt of the 
GW spectrum milder for the low frequencies than that 
in radiation dominated background, as we see in figure [6] 
(b). 

In the top hat approximation, there is an additional 
parameter x,., which can affect the shape of the GW 
spectrum. In figure [51 we show how the GW spectrum 
depends on the value of Xc- Recall that the parameter 
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FIG. 6: The spectrum of GWs produced in (a) radiation dominated era and (b) matter dominated era, evaluated by using the 
totally coherent approximation (dot-dashed line), the incoherent approximation (dotted line), and the top hat approximation 
(thick solid line). The thin solid line denoted by "simulation" represents the spectrum directly obtained from numerical 
simulations. 
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FIG. 7: The plot of F{k, Ti, Tf) given by eq. (|43|) as a function of k. We choose 13 = u — 1/2 and -y — 7 for radiation dominated 
background (a), and /3 — 2/3, u = 3/2 and 7 = 9 for matter dominated background (b). We computed for the case with 
Tf = 10ry~^ and r/ — 100?)"^ for each approximation. Other parameters are fixed to be Ti = 2rj~^ , At = lri~^ and Xc = 1. 



Xc determines the time interval beyond which the un- 
equal time correlator vanishes. If Xc is large, the source 
is correlated for a long time and the amplitude of GWs 
is suppressed in the higher frequency modes due to the 
presence of interferences. This is caused by the integrand 
cos(a; — y) in eq. (j39p . which has a dominant contribu- 
tion for high frequencies and gives a factor sin(a;c). In 
particular, the higher frequency modes tend to vanish in 
the limit Xc — >■ tt. On the other hand, if Xc is small, the 
whole amplitude of GWs is suppressed since the inter- 
val of integration in the third line of eq. (15^1) becomes 
short. In this work we fixed Xc = 1 as an intermediate 
value between two extreme cases described above. We 
note that this choice may overestimate the amplitude of 
GWs compared with the result obtained directly from 
numerical simulations in the intermediate scales between 



the Hubble radius and the width of walls. 

The advantage of approximations which we have used 
is that we can predict the slope of the GW spectrum for 
the frequencies which can not be resolved in the finite box 
lattice simulations. Assuming the top hat case, the factor 
F{k, Ti, Tf) in eq. p2)) decays like k^^ for the frequencies 
larger than ~ tJ^ in radiation dominated background. 
Note that the scale Tf corresponds to the Hubble ra- 
dius (in comoving coordinate) at the time tf. Therefore, 
the frequency at which the slope of F{k,Ti,Tf) changes 
corresponds to kh{tf) given by eq. (|^ . Furthermore, 
from eq. and table [H we expect that the factor S{k) 
in eq. (|42l) is independent of k for frequencies smaller 
than ~ kh{tf), while it decays like fc"^*^^ for frequen- 
cies larger than ~ kh{tf) and like k~^-^^ for frequencies 
larger than ~ kw{tf). Then, from eq. (|42p . we see that 
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FIG. 8: The GW spectrum evaluated by using the top hat 
approximation with varying the value of Xc- The thin solid 
line represents the spectrum directly obtained from numerical 
simulations. 



the tilt of the GW spectrum becomes for frequencies 
k < kh{tf), k^ ■ k~^ ■ cx fc~°°^ for frequencies 

kh{tf) <k< k^{tf) and k^ ■ k-^ ■ fc-si^ cx k~^-^^ for fre- 
quencies k > kfju{tf ). The fc~*''°^ behavior in the interme- 
diate frequencies represents the shghtly red tilted spec- 
trum which we found from numerical simulations in sec- 
tion llll C] However, it seems that the GW spectrum eval- 
uated by using the top hat approximation in figure|6]does 
not decay like for frequencies kh{tf) < k < kw{tf). 

This might be caused by the lack of the dynamical range. 
As we see in next subsection, when these two scales kh{tf) 
and kw(tf) become separate by many order of magni- 
tudes, the spectrum obtained by using this approxima- 
tion indeed gives the behavior in the intermediate 
scales between the Hubble radius and the width of walls. 

A similar reasoning can be applied to the case with 
matter dominated background. In matter dominated 
background, assuming the top hat case, F{k,Ti,Tf) be- 
haves like k~^ for the frequencies smaller than ^ tJ^ and 
like k"^ for the frequencies larger than ^ tJ^ . Combin- 
ing this fact with the results shown in table I, we expect 
that the tilt of the GW spectrum becomes k^ ■ fc~^ c>c k for 
frequencies k < kh{tf), k^ ■ k~^ ■ fc~^-^* cx fc^O-^s foi- fre- 
quencies kh{tf) <k<ky,(tf) and fc3-fc-i-fc-6.i2 ^ fc-4.i2 
for frequencies k > k^itf). 



D. Spectrum today 

The GW spectrum which we observe today is obtained 
by considering redshift due to the expansion of the uni- 
verse. Assuming that GWs are produced in the radiation 
dominated era, the amplitude of GWs observed today is 
given by ^] 

/100\ 

ng^h^to) = 1.67 X 10-5 i—j (rigw)*, (44) 



where is the number of relativistic degrees of freedom 
at the time when GWs are produced, h is the renormal- 
ized Hubble parameter {Hq — 100/ikmsec-^Mpc-^), and 
(ilgw)* is given by eq. ([M]) (here we use the top hat ap- 
proximation). The frequency of GWs is given by 

where oq is the scale factor at the present time, go — 3.36 
is the number of relativistic degrees of freedom today. 
To = 2.725K is the temperature of the universe observed 
today, and Tj is the temperature at the time tj. In the 
second equality of eq. (H5|) . we used a{ti) = 1, which is 
assumed in numerical simulations. Recalling that ti is 
the time at which tiT] = 1 is satisfied, we obtain 

"■<^)""(to4^)"G)- 

,(46) 

The time when the production of GWs terminates is de- 
termined by the lifetime of the domain wall networks 
which decay due to the existence of the bias term (|23p . 
It is given by [H, [13] 

tf^tdec = l^{er,r\ (47) 

Combining eqs. ^9^, dM]), gH) and (gT]), we can cal- 
culate the spectrum of GWs observed today. Strictly 
speaking, it is inappropriate to put tirj — 1, because the 
occurrence of the domain wall networks would be much 
later if the phase transition is driven by the finite temper- 
ature effect at r ~ ry. However, it hardly affects the final 
shape of the GW spectrum if domain walls survived for 
sufficiently long time, since the dominant contribution for 
the spectrum of GWs mainly come from GWs produced 
at the later time, as we discussed in section UlI Al 

We show the result of the GW spectrum calculated 
by using the procedure described above in figure 13 We 
also show the expected sensitivity from future GW ob- 
servations such as Advanced LIGO [Hi, LCGT [H], 
ET^5], LISA ;.36.], and DECIGO [s^. As we anticipated 
in [17], GWs produced by domain walls with energy scale 
T] « lO^^GeV and sufficiently small e are relevant to fu- 
ture GW direct detection experiments. In particular, the 
frequency corresponding to the Hubble radius at the de- 
cay of domain walls, at which the tilt of the spectrum 
changes, is located within the range of DECIGO and 
ground-based interferometers such as advanced LIGO, 
LCGT and ET. 

We emphasize that there is a large change in dy- 
namical range between the predicted spectra in figure [5] 
(e ^ 10-^'^-10~^'') and the spectra in figure [H] obtained 
from numerical simulations (e ^ lO^'^). We have to as- 
sume that the spectrum is fiat enough to extend over a 
large frequency range in order to give observable spectra 
which we plot in figure [SI Therefore, the results shown 
in figure [5] should be regarded as just an extrapolation of 
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the flat spectrum obtained from mimerical simulations. 
This flatness property of the GW spectra must be con- 
firmed by the future numerical studies with much larger 
dynamical range. 



V. CONCLUSIONS 

In this work, we have computed the spectrum of GWs 
produced by domain walls based on the three dimensional 
lattice simulation of the scalar field. For the evaluation 
of the GW spectrum, we apply two methods: One is 
to calculate the GW spectrum directly from numerical 
simulations as in section IIII C[ and another is to calcu- 
late it indirectly by estimating the equal time anisotropic 
stress power spectrum and using an approximation for 
the unequal time correlator as in section ITVl In the later 
case, we evaluate the GW spectrum according to the 
following procedure: First, we compute the equal time 
anisotropic stress power spectrum as eq. (P7)) . Then we 
assume the form of the power spectrum as eq. (j34p . and 
determine the fc-dependence of it by fitting with the re- 
sult of the numerical computation as eq. psp. For the 
unequal time correlator, we apply three approximations 
given by eqs. Finally we obtain the GW spec- 

trum by using the formula ((20)) . We find that the result 
with the top hat approximation given by eq. (j30[) agrees 
with the spectrum directly obtained from lattice simula- 
tions within a factor of C(l). 

From these analysises, we find the following features 
about the spectrum of GWs produced by domain walls: 

(1) The slope of the spectrum changes at two character- 
istic frequencies corresponding to the Hubble radius 
at the decay of domain walls and the width of domain 
walls. 

(2) The spectrum between two characteristic frequencies 
described above becomes flat (or slightly red tilted) 
as fc"*''^^ (in radiation dominated background). 

(3) Around the time when domain walls collapse, high 
frequency modes grow since the false vacuum regions 
fragment into small pieces (this effect is not included 
in the analysis performed in section ITVT) . 



The indirect calculation performed in section IIVI en- 
ables us to evaluate the GW spectrum beyond the 
frequencies accessible in numerical simulations. This 
method might be more effective than the direct numerical 
calculation for the problem which require a long dynam- 
ical range, like domain walls. However, we note that the 
results obtained in section IIVI rely on various nontriv- 
ial assumptions: First, the form of the function S{k, r) 
given by eq. ([55]) should be regarded as a tentative (i.e. 
other functions may also reproduce the form of H(fc, r, r) 
obtained from numerical simulations). Second, although 
the top hat ansatz (j30p reproduces the spectrum obtained 
by numerical simulations, there is no rigorous proof for 



the validity of this approximation. Finally, in the calcu- 
lation of the GW spectrum, we always assume that the 
source suddenly appears at the time U and suddenly dis- 
appears at the time tj. In reality, the the source term 
T^'^ would gradually raise and continuously decay. This 
time dependence might affect the shape of the spectrum 
of GWs, as discussed in i24|]. Nevertheless, if the future 
numerical studies enable us to have an advanced under- 
standing about the properties of the source, this indirect 
method can be the alternative to numerical simulations 
for the evaluation of the GW spectrum. Therefore, this 
work should be viewed as a first step toward the under- 
standing of the spectrum of GWs produced by domain 
walls. 
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Appendix A: Formulation of the lattice simulations 

The lattice formulation is similar to that used in our 
previous study | l7} except the points which we described 
in section IIII Al In the numerical studies, we normal- 
ize the dimensionful quantities in the unit of rj. For 
example, (p — (p/rj, t — > trj^ etc. With this normaliza- 
tion, we solve the equation of motion for the scalar fleld 
given by eq. (pij) in the three dimensional lattice by us- 
ing the fourth order Runge-Kutta method. We put the 
periodic boundary condition in the configuration of the 
scalar fields. 

The simulations are performed in the comoving box 
with size b (in the unit of r?'^). The lattice spacing is 
5x = b/N, where N is the number of grid points (here 
we take N = 256). We choose the initial time of the 
simulation so that ti = 1 and a{ti) = 1. Then, the ratio of 
the Hubble radius to the physical lattice spacing Sxphys — 
a(t)Sx is 



H- 



Sx 



phys 



N 
b/3 



(Al) 



and the ratio of the wall width du 
physical lattice spacing is 



N 



5x 



phys 



(Ai/^T])-! to the 
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FIG. 9: The GW spectrum from domain walls for the case with [r], e) = (lO^^GeV.lO"") (dotted line), e) = (10^°GeV,10~^'') 
(dashed line), and (ry, e) = (10^^GeV,10~^"^) (dash-dotted line). Solid lines represent the rough sensitivity of planned detectors. 
Other parameters are chosen so that A = 0.1 and g, = 100. 



where /? is defined by eq. We take the final time 

and the size of the comoving box of the simulation as 
tj — 151, h = 50 for the case with the radiation domi- 
nated background and t/ = 65, 6 = 24 for the case with 
the matter dominated background. We take the value of 
coupling parameter as A = 0.1. With these values of the 
parameters, at the end of the simulation, the ratios given 
by eqs. (|AT|) and (|A2| become H^^jSx^^ys ~ 125.8 < TV, 
(5^/(5 a;phys — 1.32 for the case with radiation dominated 
background, and i/~^/(5a;phys — 64.3 < N ^ (5u,/<5a;phys — 
2.09 for the case with matter dominated background. 
Therefore, these length scales are marginally resolvable 
even at the final time of the simulation. 

We give the initial conditions so that the scalar field 
has quantum fluctuation at the initial time with correla- 
tion function in the momentum space given by 

(0(k)0(k')) = ^(27r)35(3)(k + k'), 

(<^(k)<^(k')) = ^(27r)3^(3)(k + k'). (A3) 

We used the massless fluctuations as the initial condi- 
tions, since the scalar field is near the top of potential 
barrier between two vacua at the initial time. We also 
put the momentum cutoff fccut above which all fiuctua- 
tions are set to zero in order to eliminate the unphysical 
noise which comes from high frequency modes in the field 
distributions. Here we set fccut = 1- We generate initial 
conditions in momentum space as Gaussian random am- 
plitudes satisfying eq. (IA3|) . then Fourier transform them 



into the configuration space to give the spatial distribu- 
tion of the field. When we adopt these initial conditions, 
the field distribution in the momentum space is domi- 
nated by the modes around A; ~ 1 (or fc ^ ry, since we 
normalize all dimensionful quantities in the unit of r\). 
This means that, in the real space, the field value varies 
with a characteristic length scale L ^ r;"^. This length 
scale is comparable to the Hubble radius at the initial 
time, since we take ti = \ (in the unit of r\~^\ There- 
fore, we expect that these initial conditions are likely to 
lead scaling domain wall configurations which satisfy the 
property L ^ t. However, we emphasize that these initial 
conditions are chosen just for the convenience of the nu- 
merical study. We use them to remove some difficulties 
that we find when we use the thermal initial conditions 
(see section [HI Ap . Since it seems that the scaling prop- 
erty is not so much affected by the initial field configu- 
ration and we are interested in the evolution of the field 
after the formation of scaling domain wall networks, we 
expect that the results were qualitatively unchanged if 
we used the different initial conditions. 



For the calculation of the area density of domain walls 
shown in figure [I] we use the same algorithm as we used 
in m- 
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Appendix B: Direct calculation of gravitational 

waves 

In this appendix, we describe the method which is used 
in section IIII CI to calculate the spectrum of GWs from 
lattice simulations. 

Instead of using the expression ([TTl) . we replace the 
ensemble average in eq. (jl6|) by an average over a volume 
V of the comoving box, 



Pgw 



1 



1 



327rGa4 V J (27r) 



hUr,k)h'*JT,k). (Bl) 



\3 y 



Substituting the solution hij given by eqs. © and ([TU]). 
and ignoring the terms with higher order in aH , we ob- 
tain 



Pgw 



2n^G 




d^k 1 
(27r)3 



dx'V^a{x')N^ix')T^.^{T' ,k) 



dx Vx'a{x')Jy{x)T'^j^ {t , k) 



(B2) 



where we used the approximation given by eq. ([T^ . and 
averaged over a period of the oscillation of hij (r, k) with 
time. The fraction of the energy density of GWs at the 
time t given by eq. (jl9p becomes 



^gw (^) 



2G^k 



dflh 



dx'V^a{x')N,{x')T^j'^iT', k) 



dx'V^aix')Mx')T^^'^iT',k) 



(B3) 



where fJ/j is a unit vector representing the direction of k 
and dilk = dcosOdcj). The TT part of the stress-energy 
tensor is computed by applying the projection operator 
in the momentum space 



7^;jT(r,k)=A,,-fcKfc)T,,(r,k) 

= A,y,fe;(fc){9fc'/>9;(/)}(T, k), 
1 



(B4) 



A.j,fc;(fc) = P^k{k)PJl{k) - -P,,{k)Pki{k), (B5) 



(B6) 



where k — k/|k|, and {dk4>di(j)}{T, k) is the Fourier trans- 
form of dk4>{T, x.)di4>{T, x). 



Using the formulae described above, we can compute 
the GW spectrum as follows: First, we obtain the time 
evolution of the scalar field (j>{t, x) from the lattice simu- 
lation, then we compute the TT projected stress-energy 
tensor as eqs. (|B4[) - (|B6[) . Finally we perform the time in- 
tegration in eq. (|B3p to obtain the spectrum. Note that, 
as we mentioned in section lU ^gw{t) given by eq. (|B3P 
does not represent the value which would be observed 
today. In the numerical study, we compute the quan- 
tity Qg„{tf), which represents the spectrum of GWs just 
after the production of them. In order to evaluate the 
spectrum observed today, one has to multiply the dilu- 
tion factor which caused by the expansion of the universe 
(see section [TV pp . 



Appendix C: Curve fitting 



In section IIVB| we fit the equal time anisotropic 
stress power spectrum Ii{k, Tf, Tf) obtained from numer- 
ical simulations to the expression 



Sik) 



where Kfi{Tf) and K^irf) are given by eq. (15^ . The pa- 
rameters {A, B, C, D, E) are chosen so that the quantity 




E 



k ^lnn(A;) 



-{Inn(fc) - InS'(fc)}^ 



(C2) 



is minimized (hereafter we omit the time argument of 
n(A:,r, r), assuming that the fitting is applied for the 
result with r = r/). Ylk C2l) means to sum over 

all discrete values of k. crinn(fc) is the standard deviation 
of Inn(fc), which is given by 



'^\an{k) 



q'n(fc) 
W) 



(C3) 



Noting that Ii{k) is computed as an average over the di- 
rection of k [see eq. (|27p ]. we estimate o'n(fe) as a standard 
deviation of n(k) which is computed at the lattice point 
on the shell with Ikl = k. 



'n(fe) 



1 



Y^{n{k)~u{k)y 



(C4) 



k|=fc 



where is the number of lattice point (in momentum 
space) at which |k| = fc is satisfied. 

We calculate given by eq. (|C2p with varying the 
parameters (A, i?, C, D, E) and find the optimal values 
which minimize the value of x^. We also determine the 
errors of (A, B, C, D, E) from a set of values at which 
deviates by 1 from its minimum. 
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